import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse
import scipy.sparse.linalg
from time import process_time# Heat transport - implicit time discretization
Heat transport - implicit time discretization
(supplementary lecture material for the section on the heat equation https://mbd_lectures.pages.rwth-aachen.de/cmm/)
Our goal is to solve the 1d heat equation
\[\partial_t T(z,t) = \alpha \, \partial^2_z T(z,t)\]
on the domain
\[[0,z_\mathrm{max}] \times [0,t_\mathrm{end}]\]
subject to initial conditions
\[T(z,0) = T_\mathrm{ini} = \mathrm{const}\]
and boundary conditions
\[T(0,t) = T_\mathrm{surface} \quad \text{and} \quad T(z_\mathrm{max},t) = T_\mathrm{depth}.\]
See heat_equation_explcit for a more extension introduction. When applying the explicit time integration, we can observe numerical instabilities when the time step is chosen too large.
To overcome this issue, we implement an implicit time discretization.
Our computational grid is equidistant in both space and time grid \(\{(z_j, t_i)\}\) with
\[z_j \in [0,z_\mathrm{max}], \; 0\leq j \leq M, \quad \text{and} \quad t_i \in [0,t_\mathrm{end}], \; 0\leq i \leq N \]
Resulting grid size is \(\Delta z = z_\mathrm{max}/M\) and time step \(\Delta t\). Ultimately, we are interested in how the vertical temperature profile evolves with time, hence we are interested in
\[T_j^i := T(z_j,t_i).\]
Implicit refers to the fact that the finite difference quotient approximating the second order spatial derivative is written with respect to \(T_j^{i+1}\) instead of \(T_j^i\):
\[ \frac{T_j^{i+1}-T_j^i}{\Delta t} = \alpha \frac{T_{j-1}^{i+1}-2T_j^{i+1}+T_{j+1}^{i+1}}{\Delta z^2} \]
Re-arranging yields
\[ T_j^i = -\frac{\alpha \Delta t}{\Delta z^2} \left(T_{j-1}^{i+1}-2T_j^{i+1}+T_{j+1}^{i+1}\right)+T_j^{i+1}, \]
The update step in this algorithm hence constitutes a linear system solve in each step
\[ \mathbf A \, \mathbf{T}^{i+1} = \mathbf b = \mathbf{T}^i, \]
with system matrix \(\mathbf A\) und right hand side \(\mathbf b\).
Let’s observe the resulting system matrix more closely:
\[ T_j^i = -\frac{\alpha \Delta t}{\Delta z^2} \left(T_{j-1}^{i+1}-2T_{j}^{i+1}+T_{j+1}^{i+1}\right)+T_{j}^{i+1} = -F\,\color{red}{T_{j-1}^{i+1}} + (2F+1) \color{red}{T_{j}^{i+1}}-F\, \color{red}{T_{j+1}^{i+1}} \]
Expressed in matrix notation (without boundary conditions) for \(M = 4\) this reads:
\[ \begin{pmatrix} A_{jj=00} & A_{j=0,j+1} & 0 & 0 & 0\\ \color{blue}{A_{j-1,j=1}} & \color{blue}{A_{jj=11}} & \color{blue}{A_{j=1,j+1}} & 0 & 0\\ 0 & A_{j-1,j=2} & A_{jj=22} & A_{j=2,j+1} & 0\\ 0 & 0 & A_{j-1,j=3} & A_{jj=33} & 0\\ 0 & 0 & 0 & A_{j-1,j=4} & A_{MM=44}\\ \end{pmatrix} \begin{pmatrix} \color{red}{T_{0}^{i+1} \\ T_{1}^{i+1} \\ T_{2}^{i+1}} \\ T_{3}^{i+1} \\ T_{4}^{i+1} \end{pmatrix} = \begin{pmatrix} T_0^{i} \\ \color{green}{T_{1}^{i}} \\ T_{2}^{i} \\ T_3^{i} \\ T_4^{i} \end{pmatrix} \]
For \(j = 1\) the equation, for instance, reads:
\[ \color{blue}{A_{01}} \color{red}{T_{0}^{i+1}} + \color{blue}{A_{11}} \color{red}{T_{1}^{i+1}} + \color{blue}{A_{12}} \color{red}{T_{2}^{i+1}} = \color{green}{T_1^{i}}, \]
meaning that the matrix coefficients are given by
$ A_{j-1,j} = A_{j,j+1} = -F A_{jj} = 2F+1 $
The function euler_impl(…) shows, how this can be implemented.
def euler_impl(Tini, Tsurf, Tdepth, dt, tend, M, zmax, alpha, run_time_plotting = False, version='sparse'):
"""
Tini : constant initial condition
Tsurf : surface temperature
Tdepth : temperature at lower boundary
dt : constant time step size
tend : end time
M : number of spatial grid points (M+1)
zmax : extent of the computational domain
alpha : thermal diffusivity
run_time_plotting
version : dense or sparse
"""
# measuring CPU time
t0 = process_time()
# determine spatial grid resolution and initialize spatial grid
dz = zmax/M
z = np.linspace(0, zmax, M+1)
# determine number of time steps and initialize time grid
N = int(round(tend/float(dt)))
t = np.linspace(0, N*dt, N+1)
# make sure dz and dt are compatible with z and t
dz = z[1] - z[0]
dt = t[1] - t[0]
print('Spatial grid resolution: ', dz)
print('Time step: ', dt)
# determine mesh Fourier number
F = alpha * dt / dz**2
print('Mesh Fourier number: ', round(F,5))
# instantiate temperature vector
T = np.zeros(M+1)
T_i = np.zeros(M+1)
# set initial condition
for j in range(0,M+1):
T_i[j] = Tini
# set up data structure for the linear system
if version == 'dense':
# initalize dense matrix A and right hand side b
A = np.zeros((M+1, M+1))
b = np.zeros(M+1)
# define dense matrix A
for j in range(1, M):
A[j,j-1] = -F # lower # <--- here: code was to be completed
A[j,j+1] = -F # upper
A[j,j] = 1. + 2*F # main
# Dirichlet boundary condition at z = 0
A[0,0] = 1.
A[0,1] = 0
# Dirichlet BC at z = zmax
A[M,M] = 1.
A[M,M-1] = 0
# check A
# print(A)
elif version == 'sparse':
# initalize sparse matrix A diagonals and right hand side b
main = np.zeros(M+1)
lower = np.zeros(M)
upper = np.zeros(M)
b = np.zeros(M+1)
# define diagonal entries
main[:] = 1 + 2*F # <--- here: code was to be completed
lower[:] = -F
upper[:] = -F
# Dirichlet BC at z = 0
main[0] = 1.
upper[0] = 0.
# Dirichlet BC at z = zmax
main[M] = 1.
lower[M-1] = 0
A = scipy.sparse.diags(diagonals=[main, lower, upper], \
offsets=[0, -1, 1], shape=(M+1, M+1),format='csr')
# check A
# print(A)
# print(A.todense())
else:
raise ValueError('version=%s' % version)
# check if we want to vizualize during run time
if run_time_plotting is True:
fig, ax = plt.subplots()
ax.set(xlabel='T(z) [°C]', ylabel='z [m]')
ax.grid()
ax.plot(T_i,z, label='temperature profile at t = 0 s')
fig.canvas.draw()
plt.pause(1.)
# loop over all time steps
for i in range(0, N):
# compute b and solve linear system
if version == 'dense':
b = T_i
b[0] = Tsurf
b[M] = Tdepth
T[:] = scipy.linalg.solve(A,b) # <--- here: code was to be completed
elif version == 'sparse':
b = T_i
b[0] = Tsurf
b[M] = Tdepth
T[:] = scipy.sparse.linalg.spsolve(A,b) # <--- here: code was to be completed
else:
raise ValueError('version=%s' % version)
# check if we want to vizualize during run time
if run_time_plotting is True:
plt.gca().cla()
plt.gca().invert_yaxis()
ax.set(xlabel = 'T(z) [°C]', ylabel = 'z [m]')
ax.grid()
ax.plot(T, z, label = 'temperature profile at t = ' + str(dt*(i+1)) + ' s')
ax.legend()
fig.canvas.draw()
plt.pause(.1)
# switch variables before next step
T_i = T
t1 = process_time()
# check if we want to vizualize during run time
if run_time_plotting is False:
fig, ax = plt.subplots()
plt.gca().invert_yaxis()
ax.set(xlabel = 'T(z) [°C]', ylabel = 'z [m]')
ax.grid()
ax.plot(T, z, label = 'temperature profile at t = ' + str(dt*N) + ' s')
ax.legend()
fig.canvas.draw()
t1 = process_time()
print('Process time [s]', round(t1-t0,3))
return z,Talpha = 1.203e-6 # thermal diffusivity for ice in m²/s (source: https://de.wikipedia.org/wiki/Temperaturleitfähigkeit)
Tini = -5.
Tsurf = -10.
Tdepth = 0.
zmax = .1
M = 30 # 1000
dt = 5.
tend = 750.
run_time_plotting = True
version = 'sparse'
%matplotlib notebook
euler_impl(Tini, Tsurf, Tdepth, dt, tend, M, zmax, alpha, run_time_plotting, version)Spatial grid resolution: 0.0033333333333333335
Time step: 5.0
Mesh Fourier number: 0.54135
Process time [s] 11.219
(array([0. , 0.00333333, 0.00666667, 0.01 , 0.01333333,
0.01666667, 0.02 , 0.02333333, 0.02666667, 0.03 ,
0.03333333, 0.03666667, 0.04 , 0.04333333, 0.04666667,
0.05 , 0.05333333, 0.05666667, 0.06 , 0.06333333,
0.06666667, 0.07 , 0.07333333, 0.07666667, 0.08 ,
0.08333333, 0.08666667, 0.09 , 0.09333333, 0.09666667,
0.1 ]),
array([-10. , -9.64690497, -9.29467375, -8.94413236,
-8.59603294, -8.25102085, -7.90960647, -7.5721429 ,
-7.23881053, -6.9096092 , -6.58435821, -6.26270423,
-5.94413678, -5.62801054, -5.31357353, -5. ,
-4.68642647, -4.37198946, -4.05586322, -3.73729577,
-3.41564179, -3.0903908 , -2.76118947, -2.4278571 ,
-2.09039353, -1.74897915, -1.40396706, -1.05586764,
-0.70532625, -0.35309503, 0. ]))